Quantum chaos: An introduction via chains of interacting spins 1/2 
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We introduce aspects of quantum chaos by analyzing the eigenvalues and the eigenstates of quan- 
tum many-body systems. The properties of quantum systems whose classical counterparts are 
chaotic differ from those whose classical counterparts are not chaotic. The spectrum of the first 
exhibits repulsion of the energy levels. This is one of the main signatures of quantum chaos. We 
show how level repulsion develops in one-dimensional systems of interacting spins 1/2 which are 
devoid of random elements and involve only two-body interactions. In addition to the statistics of 
the eigenvalues, we analyze how the structure of the eigenstates may indicate chaos. The programs 
used to obtain the data are available online. 



I. INTRODUCTION 

Classical chaos is related to the extreme sensitivity of 
the dynamics of a system to its initial conditions, a con- 
cept that can be traced back to Poincare. 1 The main fea- 
tures of classical chaos can be illustrated by a dynamical 
billiard, which is an idealized billiard table with no fric- 
tion where a particle reflects elastically from boundaries 
which can have any shape. The motion of the particle 
is represented in phase space by a trajectory whose evo- 
lution is restricted to a surface of constant energy. De- 
pending on the shape of the boundaries, the system may 
be chaotic, which means that two trajectories whose ini- 
tial conditions are very close will diverge exponentially 
in time. The rate of this separation is characterized by 
the Lyapunov exponent.— The trajectories may also be- 
come ergodic, which implies that after a long time the 
particle will have visited the entire surface of constant 
energy. Equivalently, we may say that after a long time, 
the particle is equally likely to be found in any point of 
the accessible phase space. 

For quantum systems the notion of phase-space trajec- 
tories loses its meaning as can be seen from the Heisen- 
berg uncertainty principle. Nevertheless, because classi- 
cal physics is a limit of quantum physics, it is natural 
to search for quantum signatures of classical chaos. The 
relation between quantum mechanics and classical chaos 
is the subject of quantum chaos 

Quantum chaos is a very broad field. Our intention is 
not to review the entire field, but to focus on a particular 
aspect of it, namely, the spectral statistics and structures 
of eigenstates of quantum many-body systems as they 
transition from the regular to the chaotic domain. As a 
preparation for our goal, we first give a brief idea of some 
of the important topics of the field (more details are in 
Refs. [6H9| and in the references to be cited). 

Quantum chaos encompasses semiclassical theories 
aiming at establishing direct links between the trajec- 
tories of classically chaotic systems and the properties 
of these systems in the quantum domain. The main ad- 
vances in this area were achieved by periodic-orbit theory, 
which provides a way to calculate the spectrum of a quan- 
tum system from its classical periodic orbits (trajectories 
that repeat themselves after a certain time). An excel- 



lent introduction to the subject was given by Gutzwiller, 4 
Gutzwiller's method allowed for the development of a 
theory of "scars, "iLIS which refer to the structure of eigen- 
states that concentrate along the classical periodic orbits 
of chaotic systems. Scars were systematically studied by 
Heller jii and the first experimental observations occurred 
in quantum billiards. 3 Quantum billiards obey the laws 
of quantum mechanics and correspond to miniature ver- 
sions of dynamical billiards. 

Another important step in the understanding of quan- 
tum chaos came from the verification that the distribu- 
tion of the spacings between neighboring energy levels 
of a quantum billiard depend on the billiard's classical 



counterpart 
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ff the latter is chaotic, the energy lev- 



els are highly correlated and repel each other; if it is 
regular (integrable) , the energy levels are uncorrelated, 
randomly distributed, and can cross. 

Level repulsion is one of the main features of quantum 
chaos. It has been observed in other quantum systems, 
such as atoms in strong magnetic fields^ and systems of 
coupled particles^— Interestingly, in the context of 
number theory, level repulsion has been associated also 
with distributions of prime numbers*^ 

In this article we focus on quantum many-body sys- 
tems. A common approach when dealing with such sys- 
tems is to ignore the details of the interactions and treat 
them statistically with random matrices. The idea is 
that when the interactions are strong and the behavior 
of the system is sufficiently complex, generic properties 
should emerge. This approach was taken by Wigner— to 
describe the spectrum of heavy nuclei. He employed ma- 
trices with random elements whose only constraint was 
to satisfy the symmetries of the system. The level spac- 
ing distributions of these matrices showed level repulsion 
and agreed surprisingly well with the data from actual 
nuclei spectra. When level repulsion was later verified 
in billiards, a connection between quantum chaos and 
random matrices became established. Soon after the in- 
troduction of random matrices in nuclear physics, they 
were employed in the analysis of the spectrum of other 
quantum many-body systems, such as atoms, molecules, 
and quantum dots>£~— i 15 i 18 i 19 

The application of random matrix theory is not re- 
stricted to the statistics of eigenvalues, but accommo- 
dates also studies of eigenstates. Eigenstates of random 
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matrices are pseudo-random vectors; that is, their ampli- 
tudes are random variablesj 20 ' 21 All the eigenstates are 
statistically similar, they spread through all basis vectors 
with no preferences and are therefore ergodic. 

Despite the success of random matrix theory in de- 
scribing spectral statistical properties, it cannot capture 
the details of real quantum many-body systems. The fact 
that random matrices are completely filled with statisti- 
cally independent elements implies infinite-range inter- 
actions and the simultaneous interaction of many parti- 
cles. Real systems have few-body (most commonly only 
two-body) interactions which are usually finite range. A 
better picture of systems with finite-range interactions 
is provided by banded random matrices, which were also 
studied by Wigner. 2 - 2 - Their off-diagonal elements are ran- 
dom and statistically independent, but are non-vanishing 
only up to a fixed distance from the diagonal. There are 
also ensembles of random matrices that take into account 
the restriction to few body interactions, so that only the 
elements associated with those interactions are nonzero; 
an example is the two-body-random-ensemble 2 ^— (see 
reviews in Refs. [2lll26h . Other models which describe sys- 
tems with short-range and few-body interactions do not 
include random elements, such as nuclear shell models^ 
and the systems of interacting spins which we consider 
in this article. 

All the matrices we have mentioned can lead to level re- 
pulsion, but differences are observed. For instance, eigen- 
states of random matrices are completely spread (delo- 
calized) in any basis, whereas the eigenstates of systems 
with few-body interactions delocalize only in the middle 
of the spectrum^ - — 

In this paper we study a one-dimensional system 
of interacting spins 1/2. The system involves only 
nearest-neighbor interactions, and in some cases, also 
ncxt-ncarcst-ncighbor interactions. Depending on the 
strength of the couplings, the system may develop chaos, 
which is identified by calculating the level spacing distri- 
bution. We also compare the level of delocalization of the 
eigenstates in the integrable and chaotic domains. It is 
significantly larger in the latter case, where the most de- 
localized states are found in the middle of the spectrum. 

The paper is organized as follows. Section II provides 
a detailed description of the Hamiltonian of a spin 1/2 
chain. Section III explains how to compute the level 
spacing distribution and how to quantify the level of de- 
localization of the eigenstates. Section IV shows how 
the mixing of symmetries may erase level repulsion even 
when the system is chaotic. Final remarks are given in 
Sec. V. 

II. SPIN-1/2 CHAIN 

We study a one-dimensional spin 1/2 system (a spin 
1/2 chain) described by the Hamiltonian 

H = H Z + H NN , (la) 
where 



H, = "iSf = £ "Sf + e d S* d (lb) 

i=l \i=l / 

L-l 

H NN = ]T [J xy (S?S* +1 + SfSf +1 ) + J Z S?S! +1 ] . (lc) 

i=l 

We have set K equal to 1, L is the number of sites, 
gx, v ,z _ ^x,y,z 12 are S pj n operators at site i, and 

a*' v ' z are the Pauli matrices. The term H z gives the 
Zeeman splitting of each spin i, as determined by a static 
magnetic field in the z direction. All sites are assumed 
to have the same energy splitting u>, except a single site 
d, whose energy splitting cu + td is caused by a magnetic 
field slightly larger than the field applied on the other 
sites. This site is referred to as a defect. 

A spin in the positive z direction (up) is indicated by 
| f) or by the vector (*) ; a spin in the negative z direction 
(down) is represented by | -J.) or (J) . An up spin on site 
i has energy +Wi/2, and a down spin has energy —uji/2. 
A spin up corresponds to an excitation. 

The second term, i?NN, is known as the XXZ Hamilto- 
nian. It describes the couplings between nearest-neighbor 
(NN) spins; J xy is the strength of the flip-flop term 
SfSf +1 + SfSf +1 , and J z is the strength of the Ising inter- 
action SfSf<i- The flip-flop term exchanges the position 
of neighboring up and down spins according to 

J*y{S?Sf +1 + SfSf +1 )\ Ui+i) = {J*y/2)\ Uti+i), (2) 

or, equivalcntly, it moves the excitations through the 
chain. We have assumed open boundary conditions as 
indicated by the sum in i?NN which goes from i = 1 to 
L — X. Hence, an excitation in site 1 (or L) can move only 
to site 2 (or to site L — X). Closed boundary conditions, 
where an excitation in site 1 can move also to site L (and 
vice-versa) are mentioned briefly in Sec. IV. 

The Ising interaction implies that pairs of parallel spins 
have higher energy than pairs of anti-parallel spins, that 
is, 

■WS'+il Ui+i) = +(J,/4)| titi+i), (3) 

and 

JzS?S? +1 \ U+i) = -(J,/4)| U+i). (4) 

For the chain described by Eqs. |T]) the total spin in 
the z direction, S z = Yli=i Sf> i s conserved, that is, 
[H, S z ] = 0. This condition means that the total number 
of excitations is fixed; the Hamiltonian cannot create or 
annihilate excitations, it can only move them through the 
chain. 

To write the Hamiltonian in matrix form and diagonal- 
ize it to find its eigenvalues and eigenstates, we need to 
choose a basis. The natural choice corresponds to arrays 
of up and down spins in the z direction, as in Eqs. ©, 
© and We refer to it as the site basis. In this basis, 
H z and the Ising interaction contribute to the diagonal 
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elements of the matrix, and the flip-flop term leads to 
the off-diagonal elements. 

In the absence of the Ising interaction, the excitations 
move freely through the chain. In this case the eigenval- 
ues and eigenstates can be found analytically. The exis- 
tence of an analytical method to find the spectrum of a 
system guarantees its integrability. The addition of the 
Ising interaction eventually leads to the onset of quan- 
tum chaos. The source of chaos is the interplay between 
the Ising interaction and the defect^ 

To bring the system to the chaotic regime, we set 
Jxy — 1 (arbitrary units), choose J z = td = 0.5 (ar- 
bitrary units), and place the defect on site d = [L/2\, 
where [L/2\ stands for the largest integer less than or 
equal to L/2. These choices are based on the following 
factors, (a) The strength of the Ising interaction can- 
not be much larger than J xy , because if it were, basis 
vectors with different numbers of pairs of parallel spins 
would have very different energies and J xy would not be 
able to effectively couple them. For example, the energy 
difference between | ttt) an d | tit) is ^/ 2 ; if we had 
Jz 3> Jxy, the matrix element J xy /2 coupling these two 
basis vectors would become ineffective. As a result, the 
eigenstates would involve only a small portion of the ba- 
sis vectors, which would mean localized eigenstates and 
therefore non-chaotic systems, (b) The defect cannot be 
placed on the edges of the chain, because it has been 
shown^ that in this case an analytical solution exists 
and the system is therefore still integrable. (c) We can- 
not have ed 3> Jxy, because this would break the chain in 
two; that is, an excitation on one side of the chain would 
not have enough energy to overcome the defect and reach 
the other side of the chain. In effect, we would be deal- 
ing with two independent chains described by the XXZ 
model, which is integrable i 33 ' 34 



III. QUANTUM CHAOS 

We use the level spacing distribution to identify when 
the system becomes chaotic. We analyze also what hap- 
pens to the structure of the eigenstates once level repul- 
sion occurs. 



A. Level spacing distribution 

The distribution P(s) of the spacings, s, of neighbor- 
ing energy levels differs depending on the regime of the 
system<£&i£ The energy levels of integrable systems are 
not correlated, and are not prohibited from crossing, so 
the distribution is Poissonian (P), 

P P ( S ) = e- s . (5) 

In chaotic systems the eigenvalues become correlated 
and crossings are avoided. There is level repulsion, and 
P(s) is given by the Wigner-Dyson (WD) distribution, 



as predicted by random matrix theory. The form of 
the Wigner-Dyson distribution depends on the symme- 
try properties of the Hamiltonian. Systems with time 
reversal invariance are described by a Gaussian orthogo- 
nal ensemble, which corresponds to an ensemble of real 
symmetric matrices, whose elements are independent 
random numbers chosen from a Gaussian distribution. 
The average of the elements and the variance satisfy 
(Hij) = and (Hfj) = 1 + The level spacing dis- 
tribution of a Gaussian orthogonal ensemble is given by 

iW*) = ye-<" 2/4 - (6) 

The Hamiltonian in Eqs. ([TJ is also real and symmet- 
ric, so the distribution achieved in the chaotic limit is the 
same as Pwd(s)- However, our system has only short- 
range-two-body interactions, so it cannot reproduce all 
features of the Gaussian orthogonal ensemble. The den- 
sity of states, for instance, is Gaussian for Eqs. (JXJ) , but 
for the Gaussian orthogonal ensemble it has the shape of 
half of a circle. Interested readers can obtain these forms 
using the codes provided.— Another difference is the 
structure of the eigenstates, as discussed in Sec. lIIIBl For 
Gaussian orthogonal ensembles, all the eigenstates are 
statistically similar and highly delocalized, whereas for 
the Hamiltonian in Eqs. ([1]), derealization is restricted 
to the middle of the spectrum. 

To obtain the level spacing distribution, we first need 
to separate the eigenvalues according to their symme- 
try sectors (subspaces). If we mix eigenvalues from dif- 
ferent symmetry sectors, we cannot achieve a Wigner- 
Dyson distribution even if the system is chaotic because 
eigenvalues from different subspaces are independent and 
therefore uncorrelated, so do not repel each other. We 
discuss further the danger of mixing eigenvalues from dif- 
ferent symmetry sectors in Sec. II VI For the moment, we 
need to remember only that the Hamiltonian in Eqs. ([T]) 
conserves S z . This symmetry implies that the Hamil- 
tonian matrix is separated into uncoupled blocks, corre- 
sponding to subspaces with a fixed number of spins in 
the up direction. In the following, we select a particular 
subspace with L/3 up spins, whose Hamiltonian matrix 
has therefore dimension D = L\/[(L/3)\(L — L/3)!]. 

The second essential step before computing the spac- 
ing distribution is to unfold the spectrum. The procedure 
consists of locally rescaling the energies, so that the lo- 
cal density of states of the renormalized eigenvalues is 
1. This rescaling allows for the comparison of spectra 
obtained for different parameters and for different sys- 
tems. There are different ways to unfold the spectrum. 
A simple and commonly used procedure is to order the 
spectrum in increasing values of energy; separate it into 
several smaller sets of eigenvalues; and divide each eigen- 
value by the mean level spacing of its particular set. The 
mean level spacing of the new set of renormalized energies 
becomes 1. Because the density of states is the number 
of states in an interval of energy, that is, the reciprocal of 
the mean level spacing, this procedure also ensures that 
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the local density of states is unity. This procedure is the 
one we used^ 

Given the unfolded spacings of neighboring levels, the 
histogram can now be computed. To compare it with the 
theoretical curves, the distribution needs to be normal- 
ized, so that its total area is equal to 1. 

Figure [1] shows the level spacing distribution when the 
defect is placed on site 1 and on site \L/2\ . The first case 
corresponds to an integrable model and the distribution 
is a Poisson; the second case is a chaotic system, so the 
distribution is Wigner-Dyson. 



P 
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FIG. 1: (Color online) Level spacing distribution for the 
Hamiltonian in Eqs. {TJ with L = 15, 5 spins up, u = 0, 
= 0.5, J xy — 1, and J z = 0.5 (arbitrary units); bin size = 
0.1. (a) Defect on site d = l;(b) defect on site d = 7. The 
dashed lines are the theoretical curves. 



B. Number of principal components 

We now investigate how the transition from a Poisson 
to a Wigner-Dyson distribution affects the structure of 
the eigenstates. In particular, we study how delocalized 
they are in both regimes. 

To determine the spreading of the eigenstates in a par- 
ticular basis, we look at their components. Consider 
an eigenstate \ipi) written in the basis vectors as 
IV^} = J2k=i c ik\£,k}- It will be localized if it has the par- 
ticipation of few basis vectors, that is, if a few |cjfc | 2 make 
significant contributions. It will be delocalized if many 
\cik | 2 participate with similar values. To quantify this cri- 
terion, we use the sum of the square of the probabilities, 
\cik\ A (the sum of the probabilities would not be a good 
choice, because normalization implies Y] k _-, \cik\ 2 = 1), 
and define the number of principal components of eigen- 
state i asZ^S- 



The number of principal components gives the number 
of basis vectors which contribute to each eigenstate. It 
is small when the state is localized and large when the 
state is delocalized. 

For Gaussian orthogonal ensembles, the eigenstates are 
random vectors, that is, the amplitudes Cik are indepen- 
dent random variables. These states are completely de- 
localized. Complete derealization does not mean, how- 
ever, that the number of principal components is equal to 



D. Because the weights \cik\ 2 fluctuate, the average over 
the ensemble gives number of principal components ~ 
D/3.2L2S 

To study the number of principal components for 
Eqs. ([T]), we need to choose a basis. This choice depends 
on the question we want to address. We consider two 
bases, the site- and mean-field basis. The site-basis is 
appropriate when analyzing the spatial derealization of 
the system. To separate regular from chaotic behavior, 
a more appropriate basis consists of the eigenstates of 
the integrable limit of the model, which is known as the 
mean-field basisi^l In our case the integrable limit corre- 
sponds to Eqs. (P) with J xy ^ 0, e<j ^ 0, and J z = 0. 

We start by writing the Hamiltonian in the site-basis. 
Let us denote these basis vectors by \<j)j). In the absence 
of the Ising interaction, the diagonalization of the Hamil- 
tonian leads to the mean-field basis vectors. They are 
given by = Y]f—i bkj\<Pj}- The diagonalization of the 
complete matrix, including the Ising interaction, gives 
the eigenstates in the site-basis, \ipi) = Ylf=i a ij \4>j)- If 
we use the relation between \<f>j) and we may also 
write the eigenstates of the total Hamiltonian in Eqs. (fTJ) 
in the mean-field basis as 

d / d \ D 

= E E a ^- i&> = E c ^>- ( 8 ) 
k=i \j=i j k=i 

Figures [2] shows the number of principal components 
for the eigenstates in the site-basis [(a), (b)] and in the 
mean- field basis [(c), (d)] for the cases where the defect 
is placed on site 1 [(a), (c)] and on site \_L/2\ [(b), (d)]. 
The level of derealization increases significantly in the 
chaotic regime. However, contrary to random matrices, 
the largest values are restricted to the middle of the spec- 
trum, the states at the edges being more localized. This 
property is a consequence of the Gaussian shape of the 
density of states of systems with two-body interactions. 
The highest concentration of states appears in the middle 
of the spectrum, where the strong mixing of states can 
occur leading to widely distributed eigenstates. 

An interesting difference between the integrable and 
chaotic regimes is the fluctuations of the number of prin- 
cipal components. For the regular system the number of 
principal components shows large fluctuations. In con- 
trast, in the chaotic regime the number of principal com- 
ponents approaches a smooth function of energy. Chaotic 
eigenstates close in energy have similar structures and 
consequently similar values of the number of principal 
components. 

IV. SYMMETRIES 

The presence of a defect breaks symmetries of the sys- 
tem. In this section we remove the defect and have a 
closer look at the symmetries. 

We refer to the system in the absence of a defect 
(ed = 0) as defect-free. Contrary to the case where 
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FIG. 2: (Color online) Number of principal components for 
the eigenstates of the Hamiltonian in Eqs. (TTJ versus energy; 
L = 18, 6 spins up, uj = 0, td = 0.5, J xy = 1, and J z = 0.5 
(arbitrary units), (a) and (b) site-basis; (c) and (d) mean- 
field basis; (a) and (c) defect on site d = 1; (b) and (d) defect 
on site d = 9. 



£<i 7^ 0, a defect-free spin-1/2 chain with NN couplings 
remains integrable even when the Ising interaction is 
added. This system can be analytically solved using the 
Bethe ansatzj 33 i 34 To drive the system to chaos, while 
keeping it defect free, we need to add further couplings. 
By considering couplings between next-nearest-neighbors 
(NNNs)r^ - — the Hamiltonian becomes 



H = H] 



where 



H 



NNN 



L-2 

E 

n=l 



NN 



QX QX 

xy \ n n+2 



NNN, 



ey qy \ 
°n J n+2) 



(9) 



(10) 

For sufficiently large a (a > 0.2 for L = 15), there are 
various scenarios for which chaos can develop, which in- 
clude the absence of Ising interactions, J z = J' z = 0; 
the absence of the flip-flop term between next-nearest- 
neighbors, J' xy = 0; the absence of Ising interaction be- 
tween next-nearest-neighbors, J' z = 0; and the presence 
of all four terms. 

Depending on the parameters in Eq. , we might not 
obtain a Wigner-Dyson distribution even if the system is 
chaotic if not all symmetries of the system are taken into 
account. 37 ' 39 We have mentioned conservation of total 
spin in the z direction. In the absence of a defect other 
symmetries of H in Eq. ([§]) include the following. 40 

Parity. Parity may be understood by imagining a mir- 
ror in one edge of the chain. For eigenstates written in 
the site-basis, the probability of each basis vector is equal 
to that of its reflection. For example, suppose we have 
L = 4 and one excitation. The eigenstates are given by 
|V0 = m\ till) + a l2 | ItU) + a i3 \ UU) + (Hi\ W4t)< 
The amplitudes are either an — an and a& — aa for 
even parity or a,i = —an and — —a^ for odd parity. 
The level spacing distribution needs to be independently 
obtained for each parity. 



Spin reversal. If the chain has an even number of sites 
and L/2 up spins, then S z = 0. In this sector pairs of 
equivalent basis vectors correspond to those which be- 
come equal if we rotate all the spins from one vector by 
180°. For example, state | tilt) pairs with state | tttt)- 

Total spin. If the system is isotropic, that is, J xy — 



the total spin, S 2 = (X) n =i ^™) 2 



is 



J 2 and J xy — ^ z 
conserved. 

For closed boundary conditions there is also momen- 
tum conservation. The more symmetries the system has, 
the smaller the subspaces become for a given system size, 
which is not good for statistics. For this reason we chose 
open boundary conditions. 
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FIG. 3: (Color online) Level spacing distribution for Eq. ([9|) 
with a = 0.5, J xv = 



up, 
S z 
/' 

,J xy 



V - 7 — V - 



bin size = 0.1. (a) L = 14, 7 spins 
1. All eigenvalues of the subspace 
are considered, (b) J' xy = 1, J z = J' z = 0.5. (c) 
1, Jz = J' z = 0. (d) J' xy = 0, J z = J' z = 0.5. (e) 
J' xy = 1, Jz = 0.5, J' z = 0. For (b)-(e) L = 15, 5 spins 
up. The eigenvalues are separated according to the parity 
of the corresponding eigenstates. P{s) is the average of the 
distributions of the two parity sectors. 

In Fig. [3] we show the level spacing distribution for 
the four chaotic systems we have described. All the fig- 
ures involve eigenvalues of a single selected S^-sector. 
Both Figs.^a) and^b) show results for a chaotic spin- 
1/2 chain with all interactions in Eq. © included. Be- 
cause Fig. [3ja) mixes the three symmetries we have dis- 
cussed, P(s) becomes a Poisson distribution, even though 
the system is in fact chaotic. In Fig. [3Jb), just as in 
Figs. [3jc)-(e), we circumvent the S z — subspace by 
avoiding chains with an even number of sites. We also 
choose J z 7^ J xy to avoid conservation of total spin. By 
doing so, the only remaining symmetry is parity, which 
we take into account. The expected Wigner-Dyson dis- 
tributions are then obtained. 



V. DISCUSSION 

The computer programs used to obtain the data for 
Figs. [U [21 and [3] are available^ The reader will also find 
programs to study Gaussian orthogonal ensembles and 
suggestions for other investigations. 
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Spin 1/2 chains are excellent models for introducing 
students to some of the basic concepts of linear algebra, 
quantum mechanics, as well as to current areas of re- 
search. In addition to the crossover from integrability to 
chaos, they can be used to introduce topics as diverse 
as the metal-insulator transition, quantum phase transi- 
tion, entanglement 41 , spintronics, and methods of quan- 
tum control 4 ^. They have been considered as models for 
quantum computers 4 ^ and magnetic compounds^ 4 -, and 
recently have been simulated in optical lattices^ — 
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